Skip to content

Support wraparound indexing in longitude #47

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 3 commits into
base: main
Choose a base branch
from
Draft

Conversation

dcherian
Copy link
Collaborator

Closes #26

Apparently we inherit this from affine , so no CRS needed. Does anyone understand this?

@dcherian dcherian force-pushed the remove-rioxarray branch 7 times, most recently from b57caba to fb36168 Compare June 30, 2025 23:37
@dcherian dcherian force-pushed the remove-rioxarray branch from fb36168 to d42f267 Compare July 1, 2025 01:09
@benbovy
Copy link
Member

benbovy commented Jul 1, 2025

I think this works with the tests that you added here because the affine transformation is linear (symmetrical) and also because how numpy works with negative indices.

I don't think this is will always work without something to handle wraparound, though.

indexed = ds.pipe(assign_index)
# We lose existing attrs when calling ``assign_index``.
assert_equal(ds.sel(lon=[220, 240]), indexed.sel(lon=[-140, -120]))
assert_equal(ds.sel(lon=220), indexed.sel(lon=-140))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For example

assert_equal(ds.sel(lon=[0, 0]), indexed.sel(lon=[0, 359.9]))

should pass with some proper wraparound but raises

IndexError: index 180 is out of bounds for axis 2 with size 180

for indexed.sel(lon=359.9)

Copy link
Collaborator Author

@dcherian dcherian Jul 1, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hmmm.. using method="nearest" by default is quite wrong, but CoordinateTransformIndex requires it

@benbovy
Copy link
Member

benbovy commented Jul 1, 2025

I would be nice to also have an option for shifting the range [0, 360] from/to [-180, 180].

@dcherian dcherian force-pushed the remove-rioxarray branch from d42f267 to f4cd46a Compare July 1, 2025 20:42
@scottyhq
Copy link
Collaborator

scottyhq commented Jul 2, 2025

I would be nice to also have an option for shifting the range [0, 360] from/to [-180, 180].

Initially I thought that is what this PR was addressing. I guess I didn't realize that the affine "automatically" handles out-of-range values (for certain cases if I understand @benbovy's comment above correctly)! I was expecting a function or keyword argument to sel that handles input values in the same way that a crs argument could be provided to automatically convert coordinates before selecting.

Apparently we inherit this from affine , so no CRS needed. Does anyone understand this?

I guess CRS is needed to check that coordinates are lon/lat and raise if the CRS is projected (where there is no wrapping)? https://pyproj4.github.io/pyproj/stable/api/crs/crs.html#pyproj.crs.CRS.is_projected

In my mind there are two scenarios:

1. Selection using an affine that generates coordinates in [0,360] & user supplies out-of-bounds longitude that needs to be wrapped into the [0,360] range

def wrap_coord(lon):
  # same as converting from [-180,180] to [0,360] range
  return lon % 360
da.sel(lon=-120) == ds.isel(lon=240)
da.sel(lon=365) == ds.isel(lon=5)

2. Selection on affine that generates coordinates in [-180,180] & user supplies longitude in [0, 360] range that needs to be converted

da.sel(lon=240) == da.sel(lon=-120)
def convert_to_negative_west(lon):
  return ((lon - 180) % 360) - 180

@benbovy
Copy link
Member

benbovy commented Jul 2, 2025

Naive question: is it common to have unprojected lat/lon raster datasets with a rectilinear / no-rotation affine transform? It is pretty common for model outputs but they don't use affine transforms.

@scottyhq
Copy link
Collaborator

scottyhq commented Jul 2, 2025

Naive question: is it common to have unprojected lat/lon raster datasets with a rectilinear / no-rotation affine transform? It is pretty common for model outputs but they don't use affine transforms.

I think so, yes. For example global digital elevation datasets are often distributed as Tiff tiles (possibly w/ different but regular longitude spacing as you move towards the poles):

gdalinfo '/vsicurl?empty_dir=yes&url=https://copernicus-dem-30m.s3.eu-central-1.amazonaws.com/Copernicus_DSM_COG_10_N39_00_W109_00_DEM/Copernicus_DSM_COG_10_N39_00_W109_00_DEM.tif' -json | jq .geoTransform 
[
  -109.000138888888884,
  0.0002777777777778,
  0.0,
  40.0001388888888911,
  0.0,
  -0.0002777777777778
]
image (https://dataspace.copernicus.eu/sites/default/files/media/files/2024-06/geo1988-copernicusdem-spe-002_producthandbook_i5.0.pdf)

@benbovy
Copy link
Member

benbovy commented Jul 2, 2025

FWIW Apache SIS has a special WraparoundTransform class.

I guess one option could be to provide for convenience a similar class that would interoperate with affine.Affine (could be exposed here or in the affine package eventually?) and then let Rasterix users or downstream libraries take care of setting it appropriately (i.e., chain multiple transform objects before passing the resulting Affine object to RasterIndex.from_transform()), instead of trying to determine here if wraparound has to be applied depending on the CRS. That would be one less reason of making RasterIndex CRS-aware :).

Base automatically changed from remove-rioxarray to main July 2, 2025 14:54
@dcherian
Copy link
Collaborator Author

dcherian commented Jul 2, 2025

Yeah, this is a lot more subtle than I initially realized, I won't be able to do it before SciPy. I am going to issue an alpha release now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Support wraparound indexing for longitude
3 participants